import math
import numpy as np
import logging
from datetime import datetime

# 常量定义
PI = 3.1415927
ltStartAge = 4  # 长期平均开始年龄

# 全局变量
daysInMonth = [31, 28, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31]  # 每月的天数

# 参数字典（从 Standard 3PG parameters 和 Standard 3PGpjs options 提取）
parameters = {
    "pFS2": {"E.globulus": 1, "P. radiata": 0.8, "Trial": 1, "unit": "-"},
    "pFS20": {"E.globulus": 0.15, "P. radiata": 0.0243, "Trial": 0.15, "unit": "-"},
    "aWs": {"E.globulus": 0.095, "P. radiata": 0.8, "Trial": 0.095, "unit": "-"},
    "nWs": {"E.globulus": 2.4, "P. radiata": 0.6, "Trial": 2.4, "unit": "-"},
    "pRx": {"E.globulus": 0.8, "P. radiata": 0.8, "Trial": 0.8, "unit": "-"},
    "pRn": {"E.globulus": 0.25, "P. radiata": 0.6, "Trial": 0.25, "unit": "-"},
    "gammaF1": {"E.globulus": 0.027, "P. radiata": 0.03, "Trial": 0.027, "unit": "1/months"},
    "gammaF0": {"E.globulus": 0.001, "P. radiata": 0.001, "Trial": 0.001, "unit": "1/months"},
    "tgammaF": {"E.globulus": 12, "P. radiata": 36, "Trial": 12, "unit": "years"},
    "gammaR": {"E.globulus": 0.015, "P. radiata": 0.015, "Trial": 0.015, "unit": "1/months"},
    "Tmin": {"E.globulus": 8.5, "P. radiata": 0, "Trial": -999, "unit": "deg C"},
    "Topt": {"E.globulus": 16, "P. radiata": 20, "Trial": 25, "unit": "deg C"},
    "Tmax": {"E.globulus": 40, "P. radiata": 40, "Trial": 999, "unit": "deg C"},
    "kF": {"E.globulus": 0, "P. radiata": 1, "Trial": 1, "unit": "days"},
    "SWconst0": {"E.globulus": 0.7, "P. radiata": 0.7, "Trial": 0.7, "unit": "-"},
    "SWpower0": {"E.globulus": 9, "P. radiata": 9, "Trial": 9, "unit": "-"},
    "fCalpha700": {"E.globulus": 1.4, "P. radiata": 1.4, "Trial": 1.4, "unit": "-"},
    "fCg700": {"E.globulus": 0.7, "P. radiata": 1.4, "Trial": 1.4, "unit": "-"},
    "m0": {"E.globulus": 0, "P. radiata": 0.6, "Trial": 0, "unit": "-"},
    "fN0": {"E.globulus": 1, "P. radiata": 1, "Trial": 1, "unit": "-"},
    "fNn": {"E.globulus": 0, "P. radiata": 1, "Trial": 1, "unit": "-"},
    "MaxAge": {"E.globulus": 50, "P. radiata": 50, "Trial": 50, "unit": "years"},
    "nAge": {"E.globulus": 4, "P. radiata": 4, "Trial": 4, "unit": "-"},
    "rAge": {"E.globulus": 0.95, "P. radiata": 0.95, "Trial": 0.95, "unit": "-"},
    "gammaN1": {"E.globulus": 0, "P. radiata": 0, "Trial": 0, "unit": "%/year"},
    "gammaN0": {"E.globulus": 0, "P. radiata": 0, "Trial": 0, "unit": "%/year"},
    "tgammaN": {"E.globulus": 0, "P. radiata": 36, "Trial": 12, "unit": "years"},
    "ngammaN": {"E.globulus": 1, "P. radiata": 1, "Trial": 1, "unit": "-"},
    "wSx1000": {"E.globulus": 300, "P. radiata": 160, "Trial": 300, "unit": "kg/tree"},
    "thinPower": {"E.globulus": 1.5, "P. radiata": 1.5, "Trial": 1.5, "unit": "-"},
    "mF": {"E.globulus": 0, "P. radiata": 0.2, "Trial": 0.2, "unit": "-"},
    "mR": {"E.globulus": 0.2, "P. radiata": 0.2, "Trial": 0.2, "unit": "-"},
    "mS": {"E.globulus": 0.2, "P. radiata": 0.2, "Trial": 0.2, "unit": "-"},
    "SLA0": {"E.globulus": 11, "P. radiata": 5, "Trial": 6, "unit": "m^2/kg"},
    "SLA1": {"E.globulus": 4, "P. radiata": 5, "Trial": 2.5, "unit": "m^2/kg"},
    "tSLA": {"E.globulus": 2.5, "P. radiata": 2, "Trial": 2.5, "unit": "years"},
    "k": {"E.globulus": 0.5, "P. radiata": 0.5, "Trial": 0.5, "unit": "-"},
    "fullCanAge": {"E.globulus": 3, "P. radiata": 3, "Trial": 3, "unit": "years"},
    "MaxIntcptn": {"E.globulus": 0.15, "P. radiata": 0.15, "Trial": 0.15, "unit": "-"},
    "LAImaxIntcptn": {"E.globulus": 0, "P. radiata": 0, "Trial": 0, "unit": "-"},
    "alphaCx": {"E.globulus": 0.06, "P. radiata": 0.06, "Trial": 0.055, "unit": "molC/molPAR"},
    "y": {"E.globulus": 0.47, "P. radiata": 0.47, "Trial": 0.47, "unit": "-"},
    "MinCond": {"E.globulus": 0, "P. radiata": 0, "Trial": 0, "unit": "m/s"},
    "MaxCond": {"E.globulus": 0.02, "P. radiata": 0.02, "Trial": 0.02, "unit": "m/s"},
    "LAIgcx": {"E.globulus": 3.33, "P. radiata": 3.33, "Trial": 3.33, "unit": "-"},
    "CoeffCond": {"E.globulus": 0.05, "P. radiata": 0.05, "Trial": 0.05, "unit": "1/mBar"},
    "BLcond": {"E.globulus": 0.2, "P. radiata": 0.2, "Trial": 0.2, "unit": "m/s"},
    "fracBB0": {"E.globulus": 0.75, "P. radiata": 0.5, "Trial": 0.15, "unit": "-"},
    "fracBB1": {"E.globulus": 0.15, "P. radiata": 0.1, "Trial": 0.15, "unit": "-"},
    "tBB": {"E.globulus": 2, "P. radiata": 1.5, "Trial": 5, "unit": "years"},
    "rhoMin": {"E.globulus": 0.45, "P. radiata": 0.48, "Trial": 0.5, "unit": "tm^3"},
    "rhoMax": {"E.globulus": 0.45, "P. radiata": 0.48, "Trial": 0.5, "unit": "tm^3"},
    "tRho": {"E.globulus": 4, "P. radiata": 4, "Trial": 4, "unit": "years"},
    "aH": {"E.globulus": 0, "P. radiata": 0, "Trial": 0, "unit": "-"},
    "nHB": {"E.globulus": 0, "P. radiata": 0, "Trial": 0, "unit": "-"},
    "nHN": {"E.globulus": 0, "P. radiata": 0, "Trial": 0, "unit": "-"},
    "aV": {"E.globulus": 0, "P. radiata": 0, "Trial": 0, "unit": "-"},
    "nVB": {"E.globulus": 0, "P. radiata": 0, "Trial": 0, "unit": "-"},
    "nVN": {"E.globulus": 0, "P. radiata": 0, "Trial": 0, "unit": "-"},
    "Qa": {"E.globulus": -90, "P. radiata": -90, "Trial": -90, "unit": "W/m^2"},
    "Qb": {"E.globulus": 0.8, "P. radiata": 0.8, "Trial": 0.8, "unit": "-"},
    "gDM_mol": {"E.globulus": 24, "P. radiata": 24, "Trial": 24, "unit": "gDM/mol"},
    "molPAR_MJ": {"E.globulus": 2.3, "P. radiata": 2.3, "Trial": 2.3, "unit": "mol/MJ"},
    "Output frequency": "no",
    "Clear output region": "yes",
    "Output data": ["StemNo", "WF", "WR", "WS", "StandVol", "LAI", "avDBH"]
}

# 气候数据（Tas metdata）
Tas_metdata = {
    "TMAX1": [20.9, 21.21, 18.91, 15.79, 12.38, 10.24, 9.88, 11.12, 12.76, 15.32, 16.57, 18.5],
    "TMIN1": [8.86, 9.98, 8.39, 6.91, 5, 3.11, 2.91, 3.01, 4.23, 5.05, 6.7, 8.03],
    "PRECIP1": [71.05, 62.51, 73.28, 108.89, 125.83, 118.06, 149.54, 126.51, 121.72, 126.8, 107.8, 102.83],
    "EVAP1": [3.95, 3.9, 2.7, 1.61, 1.08, 0.61, 0.7, 1.07, 1.64, 2.32, 2.98, 3.56],
    "RADN1": [23.47, 19.74, 13.68, 8.48, 5.43, 3.87, 4.65, 7.51, 10.85, 16.07, 19.45, 21.16],
    "RAIND1": [13.32, 11.35, 14.51, 17.28, 19.04, 17.18, 20.31, 20.12, 19.63, 18.69, 17.59, 16.63],
    "FROST1": [0 for _ in range(12)]
}

# 默认参数字典（基于 assignDefaultParameters）
default_parameters = {
    "MaxAge": 50,
    "nAge": 4,
    "rAge": 0.95,
    "fullCanAge": 3,
    "k": 0.5,
    "gammaR": 0.015,
    "SWconst0": 0.7,
    "SWpower0": 9,
    "fCalpha700": 1.4,
    "fCg700": 1.4,
    "MaxIntcptn": 0.15,
    "LAImaxIntcptn": 0,
    "MinCond": 0,
    "MaxCond": 0.02,
    "LAIgcx": 3.33,
    "BLcond": 0.2,
    "CoeffCond": 0.05,
    "y": 0.47,
    "Tmax": 40,
    "Tmin": 8.5,
    "Topt": 16,
    "kF": 1,
    "pFS2": 1,
    "pFS20": 0.15,
    "aWs": 0.095,
    "nWs": 2.4,
    "pRx": 0.8,
    "pRn": 0.25,
    "m0": 0,
    "fN0": 1,
    "fNn": 0,
    "alphaCx": 0.06,
    "poolFractn": 0,
    "gammaN1": 0,
    "gammaN0": 0,
    "tgammaN": 2,
    "ngammaN": 1,
    "wSx1000": 300,
    "thinPower": 1.5,
    "mF": 0,
    "mR": 0.2,
    "mS": 0.2,
    "gammaF1": 0.027,
    "gammaF0": 0.001,
    "tgammaF": 5,
    "SLA0": 11,
    "SLA1": 4,
    "tSLA": 2.5,
    "fracBB0": 0.75,
    "fracBB1": 0.15,
    "tBB": 2,
    "rho0": 0.45,
    "rho1": 0.45,
    "tRho": 4,
    "aH": 0,
    "nHB": 0,
    "nHN": 0,
    "aV": 0,
    "nVB": 0,
    "nVN": 0,
    "Qa": -90,
    "Qb": 0.8,
    "gDM_mol": 24,
    "molPAR_MJ": 2.3,
}

# 站点参数类
class SiteData:
    def __init__(self):
        self.Sites = ""
        self.Latitude = 0.0
        self.FR = 0.0
        self.Soil_type = ""
        self.Max_ASW = 0.0
        self.Min_ASW = 0.0
        self.Initial_stocking = 0.0
        self.alphaCx = 0.0
        self.gammaNx = 0.0

# 模拟类
class ThreePGModel:
    def __init__(self, site_data, species, climate_data, seedling_mass, month_planted, planted_year, end_age, lookup, params, output_frequency="m", silvicultural_events=None):
        self.site_data = site_data
        self.species = species
        self.climate_data = climate_data
        self.seedling_mass = seedling_mass
        self.month_planted = month_planted
        self.planted_year = planted_year
        self.end_age = end_age
        self.silvicultural_events = silvicultural_events or {}
        self.interpolate_lookups = lookup

        # self.params = params
        # logging.info(f"Initialized with params: {self.params}")  # 添加调试日志
        # self.output_frequency = output_frequency  # 新增：保存输出频率
        # # 更新 params 中依赖 site_data 的值
        # self.params["alphaCx"] = self.site_data.alphaCx
        # self.params["gammaN1"] = self.site_data.gammaNx
        # self.params["gammaN0"] = 0

        self.params = params.copy()  # 使用副本，避免修改原始字典
        #logging.info(f"Initialized with params: {self.params}")
        #print(f"Initial SLA0: {self.params['SLA0']}, SLA1: {self.params['SLA1']}")  # 添加：检查初始 SLA0 和 SLA1
        # 仅在 params 中缺失时使用 site_data 的值
        self.output_frequency = output_frequency  # 新增：保
        self.params["alphaCx"] = self.params.get("alphaCx", site_data.alphaCx)
        self.params["gammaN1"] = self.params.get("gammaN1", site_data.gammaNx)
        self.params["gammaN0"] = self.params.get("gammaN0", 0)


        self.climate = self._load_climate(climate_data)
        self._initialize()
        self.computeMetdataVbls()

        if self.params["MaxAge"] <= 0:
            raise ValueError("MaxAge must be greater than 0")
        if self.params["fullCanAge"] <= 0:
            raise ValueError("fullCanAge must be greater than 0")
        if self.params["tgammaF"] <= 0:
            raise ValueError("tgammaF must be greater than 0")
        if self.params["gDM_mol"] <= 0:
            raise ValueError("gDM_mol must be greater than 0")
        if self.params["molPAR_MJ"] <= 0:
            raise ValueError("molPAR_MJ must be greater than 0")

        # _load_parameters 方法不再需要，可以删除或保留供其他场景使用

    def _load_parameters(self, species):
        species_params = default_parameters.copy()
        for param, values in parameters.items():
            if isinstance(values, dict) and species in values:
                species_params[param] = values[species]
        return species_params

    def _load_climate(self, climate_data):
        climate = {}
        key_mapping = {
            "TMAX1": "TMAX1",
            "TMIN1": "TMIN1",
            "PRECIP1": "mRain",
            "EVAP1": "EVAP1",
            "RADN1": "mSolarRad",
            "FROST1": "mFrostDays"
        }
        for key, value in climate_data.items():
            mapped_key = key_mapping.get(key, key)
            climate[mapped_key] = np.array(value * (self.end_age + 1))
        climate["mYears"] = self.end_age + 1
        return climate

    def _initialize(self):
        self.month_map = {"January": 1, "February": 2, "March": 3, "April": 4, "May": 5, "June": 6,
                          "July": 7, "August": 8, "September": 9, "October": 10, "November": 11, "December": 12,
                          "Jan": 1, "Feb": 2,"Aug": 8,"Sep": 9,"Oct": 10,"Nov": 11,"Dec": 12
                        }
        self.InitialYear = self.planted_year
        self.InitialMonth = self.month_map[self.month_planted]
        self.PlantedYear = self.planted_year
        self.PlantedMonth = self.month_map[self.month_planted]
        self.getStandAge()
        self.NPP = 0
        self.cEpsilonStem=0
        self.Lat = self.site_data.Latitude
        self.FR = self.site_data.FR
        self.MaxASW = self.site_data.Max_ASW
        self.MinASW = self.site_data.Min_ASW
        self.StemNoi = self.site_data.Initial_stocking
        self.ASWi = 99999999
        self.soilClass = 3 if self.site_data.Soil_type in ["3", "CL"] else (
            0 if self.site_data.Soil_type == "0" else (
                1 if self.site_data.Soil_type in ["1", "S"] else (
                    2 if self.site_data.Soil_type in ["2", "SL"] else (
                        4 if self.site_data.Soil_type in ["4", "C"] else (
                                                                             -1 if self.site_data.Soil_type == "?" else 0
                        )))))
        if self.soilClass > 0:
            self.SWconst = 0.8 - 0.1 * self.soilClass
            self.SWpower = 11 - 2 * self.soilClass
        elif self.soilClass < 0:
            self.SWconst = self.params["SWconst0"]
            self.SWpower = self.params["SWpower0"]
        else:
            self.SWconst = 999
            self.SWpower = self.params["SWpower0"]

        self.params["pfsPower"] = math.log(self.params["pFS20"] / self.params["pFS2"]) / math.log(20 / 2)
        self.params["pfsConst"] = self.params["pFS2"] / (2 ** self.params["pfsPower"])
        
        self.params["fCalphax"] = self.params["fCalpha700"] / (2 - self.params["fCalpha700"])
        self.params["fCg0"] = self.params["fCg700"] / (2 * self.params["fCg700"] - 1)

        self.fractionWF = 0.5
        self.fractionWR = 0.25
        self.fractionWS = 0.25
        total_mass = self.seedling_mass * self.StemNoi / 1e6
        # print('290lines')
        # print(self.seedling_mass)
        print(self.StemNoi)
        self.WFi = (self.fractionWF / (self.fractionWF + self.fractionWR + self.fractionWS)) * total_mass
        self.WRi = (self.fractionWR / (self.fractionWF + self.fractionWR + self.fractionWS)) * total_mass
        self.WSi = (self.fractionWS / (self.fractionWF + self.fractionWR + self.fractionWS)) * total_mass
        # print('293lines')
        # print(total_mass)

        self.WLi = 0
        self.WS = self.WSi
        self.WF = self.WFi
        self.WR = self.WRi
        self.WL = self.WLi
        self.TotalW = self.WSi + self.WFi + self.WRi
        self.StemNo = self.StemNoi
        self.ASW = self.ASWi
        self.thinEventNo = 0  # 疏伐事件计数器，从 0 开始
        self.defoltnEventNo = 0  # 落叶事件计数器，从 0 开始
        self.ltStemGR = 0
        self.ltAge0 = max(self.StandAge, ltStartAge)
        self.pooledSW = 0

        self.SLA = self.expF(self.StandAge, self.params["SLA0"], self.params["SLA1"], self.params["tSLA"], 2)
        self.fracBB = self.expF(self.StandAge, self.params["fracBB0"], self.params["fracBB1"], self.params["tBB"], 1)
        self.Density = self.expF(self.StandAge, self.params["rho0"], self.params["rho1"], self.params["tRho"], 1)
        self.gammaF = self.gammaFoliage(self.StandAge)

        self.AvStemMass = self.WS * 1000 / self.StemNo
        # print('314lines')
        # print(self.WS)
        self.avDBH = (self.AvStemMass / self.params["aWs"]) ** (1 / self.params["nWs"])
        # print('316行')
        # print(self.AvStemMass)
        # print(self.avDBH)
        self.BasArea = (((self.avDBH / 200) ** 2) * PI) * self.StemNo
        self.Height = self.params["aH"] * (self.avDBH ** self.params["nHB"]) * (self.StemNo ** self.params["nHN"])
        self.LAI = self.WF * self.SLA * 0.1
        self.cLAI = self.LAI
        print(f"Initial StandAge: {self.StandAge}, SLA0: {self.params['SLA0']}, SLA1: {self.params['SLA1']}, SLA: {self.SLA}, WF: {self.WF}, LAI: {self.LAI}")  # 添加：初始值

        self.StandVol = self.WS * (1 - self.fracBB) / self.Density if self.params["aV"] == 0 else \
            self.params["aV"] * (self.avDBH ** self.params["nVB"]) * (self.StemNo ** self.params["nVN"])
        self.oldVol = self.StandVol
        self.MAI = self.StandVol / self.StandAge if self.StandAge > 0 else 0

        if self.ASWi <= self.MinASW:
            self.ASWi = self.MinASW
        elif self.ASWi >= self.MaxASW:
            self.ASWi = self.MaxASW
        self.ASW = self.ASWi

        self.noopAges = 0
        self.opAges = []
        self.opFrequency = 3
        self.monthsInStep = 1
        self.opDateNo = 1
        self.mortality = 0
        self.selfThin = 0
        self.aStemDM = 0
        self.aRADint = 0
        self.aGPP = 0
        self.aNPP = 0
        self.aEvapTransp = 0
        self.aTransp = 0
        self.aRunOff = 0
        self.aSupIrrig = 0
        self.LAIx = 0
        self.ageLAIx = 0
        self.MAIx = 0
        self.ageMAIx = 0
        self.applIrrig = 0

    def getStandAge(self):
        planted_date = datetime(self.PlantedYear, self.PlantedMonth, 1)
        initial_date = datetime(self.InitialYear, self.InitialMonth, 1)
        months_diff = (initial_date.year - planted_date.year) * 12 + (initial_date.month - planted_date.month)
        self.StandAge = months_diff / 12
        self.StartAge = int(self.StandAge)
        if self.StandAge < 0 or self.StandAge > self.end_age:
            raise ValueError(f"Invalid StandAge: {self.StandAge}")

    def expF(self, x, g0, gx, tg, ng):
        if tg != 0:
            # print(f"expF: x={x}, g0={g0}, gx={gx}, tg={tg}, ng={ng}")  # 添加：检查计算过程
            return gx + (g0 - gx) * math.exp(-math.log(2) * (x / tg) ** ng)
        return gx

    def gammaFoliage(self, age):
        if self.params["tgammaF"] * self.params["gammaF1"] == 0:
            return self.params["gammaF1"]
        kgammaF = 12 * math.log(1 + self.params["gammaF1"] / self.params["gammaF0"]) / self.params["tgammaF"]
        return (self.params["gammaF1"] * self.params["gammaF0"]) / \
               (self.params["gammaF0"] + (self.params["gammaF1"] - self.params["gammaF0"]) * math.exp(-kgammaF * age))

    def getConductance(self, LAI):
        if LAI <= self.params["LAIgcx"]:
            gC = self.params["MinCond"] + (self.params["MaxCond"] - self.params["MinCond"]) * LAI / self.params["LAIgcx"]
        else:
            gC = self.params["MaxCond"]
        return gC * self.PhysMod * self.fCg

    def getTranspiration(self, Q, VPD, h, gBL, gC):
        e20 = 2.2
        rhoAir = 1.2
        lambda_ = 2460000
        VPDconv = 0.000622
        netRad = self.params["Qa"] + self.params["Qb"] * (Q * 1e6 / h)
        defTerm = rhoAir * lambda_ * (VPDconv * VPD) * gBL
        div = gC * (1 + e20) + gBL
        Etransp = gC * (e20 * netRad + defTerm) / div
        return Etransp / lambda_ * h

    def getDayLength(self, Lat, dayOfYear):
        SLAt = math.sin(PI * Lat / 180)
        cLat = math.cos(PI * Lat / 180)
        sinDec = 0.4 * math.sin(0.0172 * (dayOfYear - 80))
        cosH0 = -sinDec * SLAt / (cLat * math.sqrt(1 - sinDec ** 2))
        if cosH0 > 1:
            return 0
        elif cosH0 < -1:
            return 1
        return math.acos(cosH0) / PI

    def getVPD(self, Tx, Tn):
        VPDx = 6.1078 * math.exp(17.269 * Tx / (237.3 + Tx))
        VPDn = 6.1078 * math.exp(17.269 * Tn / (237.3 + Tn))
        return (VPDx - VPDn) / 2

    def computeMetdataVbls(self):
        mDayLength = []
        mTav = []
        mVPD = []
        dayofyr = 0
        for month in range(12 * self.climate["mYears"]):
            if (month % 12) == 0:
                dayofyr = -15
            dayofyr += daysInMonth[(month % 12)]
            mDayLength.append(86400 * self.getDayLength(self.Lat, dayofyr))
            mTav.append((self.climate["TMAX1"][month] + self.climate["TMIN1"][month]) / 2)
            mVPD.append(self.getVPD(self.climate["TMAX1"][month], self.climate["TMIN1"][month]))
        self.climate["mDayLength"] = np.array(mDayLength)
        self.climate["mTav"] = np.array(mTav)
        self.climate["mVPD"] = np.array(mVPD)

    def getMortality(self, oldN, oldW):
        accuracy = 1 / 1000
        n = oldN / 1000
        x1 = 1000 * self.params["mS"] * oldW / oldN
        i = 0
        while True:
            i += 1
            x2 = self.params["wSx1000"] * (n ** (1 - self.params["thinPower"]))
            print(f"Iteration {i}: wSx1000={self.params['wSx1000']}, n={n}, x2={x2}")  # 添加：迭代中 wSx1000
            fN = x2 - x1 * n - (1 - self.params["mS"]) * oldW
            dfN = (1 - self.params["thinPower"]) * x2 / n - x1
            dN = -fN / dfN
            n += dN
            if abs(dN) <= accuracy or i >= 5:
                break
        return oldN - 1000 * n

    def lookup_table(self, age, ages, values):
        """查找表值，支持跳跃式或线性插值"""
        if not ages or not values or len(ages) != len(values):
            return None
        if age <= ages[0]:
            return values[0]
        if age >= ages[-1]:
            return values[-1]

        for i in range(len(ages) - 1):
            if ages[i] <= age < ages[i + 1]:
                if self.interpolate_lookups:
                    # 线性插值
                    t = (age - ages[i]) / (ages[i + 1] - ages[i])
                    return values[i] + t * (values[i + 1] - values[i])
                else:
                    # 跳跃式
                    return values[i]
        return values[-1]

    def doThinning(self):
        """执行疏伐事件，与 VBA 的 doThinning 逻辑一致"""
        if "Thinning" not in self.silvicultural_events or not self.silvicultural_events["Thinning"]["ages"]:
            return
        thin_data = self.silvicultural_events["Thinning"]
        ages = thin_data["ages"]
        if self.thinEventNo < len(ages) and self.StandAge >= ages[self.thinEventNo]:
            # print(f"Thinning at StandAge: {self.StandAge}")
            residual_stems = thin_data["residual_stems"][self.thinEventNo]
            thinWF = thin_data["thinWF"][self.thinEventNo]
            thinWR = thin_data["thinWR"][self.thinEventNo]
            thinWS = thin_data["thinWS"][self.thinEventNo]

            if self.StemNo > residual_stems:
                delN = (self.StemNo - residual_stems) / self.StemNo
                self.StemNo = self.StemNo * (1 - delN)
                # print(f"The value of StemNosss is: {self.StemNo}")
                self.WF *= (1 - delN * thinWF)
                self.WR *= (1 - delN * thinWR)
                self.WS *= (1 - delN * thinWS)

            self.thinEventNo += 1

    def doDefoliation(self):
        """执行落叶事件"""
        if "Defoliation" not in self.silvicultural_events or not self.silvicultural_events["Defoliation"]["ages"]:
            return
        defol_data = self.silvicultural_events["Defoliation"]
        ages = defol_data["ages"]
        if self.defoltnEventNo < len(ages) and self.StandAge >= ages[self.defoltnEventNo]:
           # print(f"Defoliation at StandAge: {self.StandAge}")
            defol_fraction = defol_data["defol_fraction"][self.defoltnEventNo]
            self.WF *= defol_fraction
            self.defoltnEventNo += 1

    def assignSilviculturalEvents(self, age):
        """分配营林事件，只更新 Fertility 和 Irrigation"""
        if "Fertility" in self.silvicultural_events:
            fert_data = self.silvicultural_events["Fertility"]
            fr = self.lookup_table(age+1, fert_data["ages"], fert_data["FR"])

            if fr is not None:
                self.FR = fr

        if "Irrigation" in self.silvicultural_events:
            irrig_data = self.silvicultural_events["Irrigation"]
            irrig = self.lookup_table(age, irrig_data["ages"], irrig_data["irrig"])
            self.applIrrig = irrig if irrig is not None else 0

    def assignVaryData(self, age):
        pass

    def run3PG(self, output_vars):
        #logging.info(f"Step 1: Using params {self.params}")
        output_mapping = {
            "StemNo": lambda self: self.StemNo,
            "WF": lambda self: self.WF,
            "WR": lambda self: self.WR,
            "WS": lambda self: self.WS,
            "WL": lambda self: self.WL,
            "StandVol": lambda self: self.StandVol,
            "LAI": lambda self: self.LAI,
            "avDBH": lambda self: self.avDBH,
            "Irrig": lambda self: self.applIrrig,
            "incrWS": lambda self: self.incrWS,
            "mortality": lambda self: self.mortality,
            "selfThin": lambda self: self.selfThin,
            "meantemp": lambda self: self.Tav,
            "FR": lambda self: self.FR,
            "CO2": lambda self: 350,
            "minASW": lambda self: self.MinASW,
            "StandAge": lambda self: self.StandAge,
            "BasArea": lambda self: self.BasArea,
            "Height": lambda self: self.Height,
            "MAI": lambda self: self.MAI,
            "MAIx": lambda self: self.MAIx,
            "ageMAIx": lambda self: self.ageMAIx,
            "SLA": lambda self: self.SLA,
            "CanCover": lambda self: self.CanCover,
            "LAIx": lambda self: self.LAIx,
            "ageLAIx": lambda self: self.ageLAIx,
            "TotalW": lambda self: self.TotalW,
            "AvStemMass": lambda self: self.AvStemMass,
            "Density": lambda self: self.Density,
            "fracBB": lambda self: self.fracBB,
            "fAge": lambda self: self.fAge,
            "fVPD": lambda self: self.fVPD,
            "fTemp": lambda self: self.fT,
            "fT": lambda self: self.fT,
            "fCalpha": lambda self: self.fCalpha,
            "fCg": lambda self: self.fCg,
            "fFrost": lambda self: self.fFrost,
            "fSW": lambda self: self.fSW,
            "fNutr": lambda self: self.fNutr,
            "PhysMod": lambda self: self.PhysMod,
            "GPP": lambda self: self.GPP,
            "NPP": lambda self: self.NPP,
            "RADint": lambda self: self.RADint,
            "alphaC": lambda self: self.alphaC,
            "Epsilon": lambda self: self.epsilon,
            "CVI": lambda self: self.CVI,
            "m": lambda self: self.m,
            "pR": lambda self: self.pR,
            "pS": lambda self: self.pS,
            "pF": lambda self: self.pF,
            "pFS": lambda self: self.pFS,
            "gammaF": lambda self: self.gammaF,
            "wSmax": lambda self: self.wSmax,
            "gammaN": lambda self: self.gammaN,
            "supIrrig": lambda self: self.SupIrrig,
            "RunOff": lambda self: self.RunOff,
            "fRainInt": lambda self: self.fracRainIntcptn,
            "RainInt": lambda self: self.RainIntcptn,
            "CanCond": lambda self: self.CanCond,
            "WUE": lambda self: self.WUE,
            "ET": lambda self: self.EvapTransp,
            "Transp": lambda self: self.Transp,
            "ASW": lambda self: self.ASW,
            "DayLength": lambda self: self.DayLength,
            "FrostDays": lambda self: self.FrostDays,
            "SolarRad": lambda self: self.SolarRad,
            "Tav": lambda self: self.Tav,
            "VPD": lambda self: self.VPD,
            "Rain": lambda self: self.Rain,
            "ltStemGR": lambda self: self.ltStemGR,
            "Litter": lambda self:self.cLitter,
            "TotalLitter": lambda self: self.WL,
            "EpsilonStem":lambda self: self.cEpsilonStem,
            "pooledSW": lambda self: self.pooledSW,
        }
        results = {var: [] for var in output_vars}
        results["Time"] = []
        metMonth = self.InitialMonth - 1
        current_year = self.InitialYear
        current_month = self.InitialMonth
        nThin = len(self.silvicultural_events.get("Thinning", {}).get("ages", []))
        nDefol = len(self.silvicultural_events.get("Defoliation", {}).get("ages", []))
        target_month = self.month_map[self.month_planted]  # 目标月份（种植月份）
        for year in range(self.end_age - self.StartAge):
            self.aStemDM = 0
            self.aRADint = 0
            self.aGPP = 0
            self.aNPP = 0
            self.aEvapTransp = 0
            self.aTransp = 0
            self.aRunOff = 0
            self.aSupIrrig = 0

            month = self.InitialMonth - 1
            for monthCounter in range(12):
                if self.opFrequency == 3 or monthCounter == 0:
                    self.selfThin = 0
                    self.mortality = 0
                    self.cRADint = 0
                    self.cLAI = 0
                    self.cCVI = 0
                    self.cGPP = 0
                    self.cNPP = 0
                    self.cStemDM = 0
                    self.cTransp = 0
                    self.cEvapTransp = 0
                    self.cRainInt = 0
                    self.cRunOff = 0
                    self.cSupIrrig = 0
                    self.cLitter = 0
                if abs(self.StandAge - self.ltAge0) < 0.001:
                    self.ltWS0 = self.WS

                self.assignSilviculturalEvents(self.StandAge)
                self.assignVaryData(self.StandAge)
                month = (month + 1) % 12
                metMonth = (metMonth + 1) % (12 * self.climate["mYears"])
                if metMonth == 0:
                    metMonth = 12 * self.climate["mYears"]

                self.SolarRad = self.climate["mSolarRad"][metMonth]
                self.Tav = self.climate["mTav"][metMonth]
                self.VPD = self.climate["mVPD"][metMonth]
                self.FrostDays = self.climate["mFrostDays"][metMonth]
                self.DayLength = self.climate["mDayLength"][metMonth]
                self.Rain = self.climate["mRain"][metMonth]

                # 计算气象和生理变量（保持不变）
                if self.Tav <= self.params["Tmin"] or self.Tav >= self.params["Tmax"]:
                    self.fT = 0
                else:
                    self.fT = ((self.Tav - self.params["Tmin"]) / (self.params["Topt"] - self.params["Tmin"])) * \
                              ((self.params["Tmax"] - self.Tav) / (self.params["Tmax"] - self.params["Topt"])) ** \
                              ((self.params["Tmax"] - self.params["Topt"]) / (
                                          self.params["Topt"] - self.params["Tmin"]))
                self.fVPD = math.exp(-self.params["CoeffCond"] * self.VPD)
                self.MoistRatio = self.ASW / self.MaxASW
                self.fSW = 1 / (1 + ((1 - self.MoistRatio) / self.SWconst) ** self.SWpower)
                self.fNutr = 1 if self.params["fNn"] == 0 else 1 - (1 - self.params["fN0"]) * (1 - self.FR) ** \
                                                               self.params["fNn"]
                self.fFrost = 1 - self.params["kF"] * (self.FrostDays / 30)
                self.fCalpha = self.params["fCalphax"] * 350 / (350 * (self.params["fCalphax"] - 1) + 350)
                self.fCg = self.params["fCg0"] / (1 + (self.params["fCg0"] - 1) * 350 / 350)
                self.fAge = 1 if self.params["nAge"] == 0 else 1 / (
                            1 + (self.StandAge / self.params["MaxAge"] / self.params["rAge"]) ** self.params["nAge"])
                self.PhysMod = min(self.fVPD, self.fSW) * self.fAge

                self.CanCover = (self.StandAge + 0.01) / self.params["fullCanAge"] if (
                            self.params["fullCanAge"] > 0 and self.StandAge < self.params["fullCanAge"]) else 1
                self.lightIntcptn = 1 - math.exp(-self.params["k"] * self.LAI / self.CanCover)
                self.alphaC = self.params["alphaCx"] * self.fNutr * self.fT * self.fFrost * self.fCalpha * self.PhysMod
                self.epsilon = self.params["gDM_mol"] * self.params["molPAR_MJ"] * self.alphaC
                self.RAD = self.SolarRad * daysInMonth[month]
                self.RADint = self.RAD * self.lightIntcptn * self.CanCover
                self.GPP = self.epsilon * self.RADint / 100
                self.NPP = self.GPP * self.params["y"]

                self.CanCond = self.getConductance(self.LAI)
                self.Transp = daysInMonth[month] * self.getTranspiration(self.SolarRad, self.VPD, self.DayLength,
                                                                         self.params["BLcond"], self.CanCond)
                self.fracRainIntcptn = self.params["MaxIntcptn"] * min(1, self.LAI / self.params["LAImaxIntcptn"]) if \
                self.params["LAImaxIntcptn"] > 0 else self.params["MaxIntcptn"]
                self.RainIntcptn = self.Rain * self.fracRainIntcptn

                self.SupIrrig = 0
                self.RunOff = 0
                self.ASW += self.Rain + (100 * self.applIrrig / 12) + self.pooledSW
                self.EvapTransp = min(self.ASW, self.Transp + self.RainIntcptn)
                self.excessSW = max(self.ASW - self.EvapTransp - self.MaxASW, 0)
                self.ASW -= self.EvapTransp + self.excessSW
                self.pooledSW = self.params["poolFractn"] * self.excessSW
                self.RunOff = (1 - self.params["poolFractn"]) * self.excessSW
                if self.ASW < self.MinASW:
                    self.SupIrrig = self.MinASW - self.ASW
                    self.ASW = self.MinASW

                self.TranspScaleFactor = self.EvapTransp / (self.Transp + self.RainIntcptn)
                self.GPP *= self.TranspScaleFactor
                self.NPP *= self.TranspScaleFactor

                self.WUE = 100 * self.NPP / self.EvapTransp if self.EvapTransp != 0 else 0
                self.m = self.params["m0"] + (1 - self.params["m0"]) * self.FR
                # print(self.avDBH)
                # 计算 pFS 时避免负幂问题
                if self.avDBH > 0 or self.params["pfsPower"] >= 0:
                    self.pFS = self.params["pfsConst"] * (self.avDBH ** self.params["pfsPower"])
                else:
                    self.pFS = self.params["pfsConst"]  # 使用默认值避免错误
                    print(f"警告：avDBH={self.avDBH}, pfsPower={self.params['pfsPower']}，设置 pFS 为 pfsConst")
                # print(self.avDBH)
                self.pR = self.params["pRx"] * self.params["pRn"] / (
                            self.params["pRn"] + (self.params["pRx"] - self.params["pRn"]) * self.PhysMod * self.m)
                self.pS = (1 - self.pR) / (1 + self.pFS)
                self.pF = 1 - self.pR - self.pS
                self.incrWF = self.NPP * self.pF
                self.incrWR = self.NPP * self.pR
                self.incrWS = self.NPP * self.pS
                self.lossWF = self.gammaF * self.WF
                self.lossWR = self.params["gammaR"] * self.WR

                self.WF += self.incrWF - self.lossWF
                self.WR += self.incrWR - self.lossWR
                self.WS += self.incrWS
                self.WL += self.lossWF
                self.TotalW = self.WF + self.WR + self.WS

                self.StandAge += 1 / 12

                if self.thinEventNo < nThin:
                    self.doThinning()
                if self.defoltnEventNo < nDefol:
                    self.doDefoliation()

                self.gammaN = self.expF(self.StandAge, self.params["gammaN0"], self.params["gammaN1"],
                                        self.params["tgammaN"], self.params["ngammaN"])
                if self.gammaN > 0:
                    delStems = self.gammaN * self.StemNo / 12 / 100
                    self.WF -= self.params["mF"] * delStems * (self.WF / self.StemNo)
                    self.WR -= self.params["mR"] * delStems * (self.WR / self.StemNo)
                    self.WS -= self.params["mS"] * delStems * (self.WS / self.StemNo)
                    self.StemNo -= delStems
                    self.mortality += delStems
                    

                self.wSmax = self.params["wSx1000"] * (1000 / self.StemNo) ** self.params["thinPower"]
                self.AvStemMass = self.WS * 1000 / self.StemNo

               
                if self.wSmax < self.AvStemMass:
                    delStems = self.getMortality(self.StemNo, self.WS)
                    self.WF -= self.params["mF"] * delStems * (self.WF / self.StemNo)
                    self.WR -= self.params["mR"] * delStems * (self.WR / self.StemNo)
                    self.WS -= self.params["mS"] * delStems * (self.WS / self.StemNo)
                    self.StemNo -= delStems
                    self.wSmax = self.params["wSx1000"] * (1000 / self.StemNo) ** self.params["thinPower"]
                    self.AvStemMass = self.WS * 1000 / self.StemNo
                    self.selfThin += delStems
                   

                self.SLA = self.expF(self.StandAge, self.params["SLA0"], self.params["SLA1"], self.params["tSLA"], 2)
                self.fracBB = self.expF(self.StandAge, self.params["fracBB0"], self.params["fracBB1"],
                                        self.params["tBB"], 1)
                self.Density = self.expF(self.StandAge, self.params["rho0"], self.params["rho1"], self.params["tRho"],
                                         1)
                self.gammaF = self.gammaFoliage(self.StandAge)

                self.LAI = self.WF * self.SLA * 0.1
                #print(f"Step {year}-{monthCounter+1}: StandAge={self.StandAge}, SLA0={self.params['SLA0']}, SLA1={self.params['SLA1']}, SLA={self.SLA}, WF={self.WF}, LAI={self.LAI}")  # 添加：动态值
                # 计算 AvStemMass 和 avDBH 时添加检查
                self.AvStemMass = self.WS * 1000 / self.StemNo if self.StemNo > 0 else 0
                if self.AvStemMass > 0 and self.params["aWs"] > 0 and self.params["nWs"] > 0:
                    self.avDBH = (self.AvStemMass / self.params["aWs"]) ** (1 / self.params["nWs"])
                else:
                    self.avDBH = 0.01  # 使用一个小的正值避免 0
                    print(
                        f"警告：AvStemMass={self.AvStemMass}, aWs={self.params['aWs']}, nWs={self.params['nWs']}，设置 avDBH 为 0.01")

                self.BasArea = (((self.avDBH / 200) ** 2) * PI) * self.StemNo
                self.Height = self.params["aH"] * (self.avDBH ** self.params["nHB"]) * (
                            self.StemNo ** self.params["nHN"])
                self.StandVol = self.WS * (1 - self.fracBB) / self.Density if self.params["aV"] == 0 else \
                    self.params["aV"] * (self.avDBH ** self.params["nVB"]) * (self.StemNo ** self.params["nVN"])

                self.CVI = self.StandVol - self.oldVol
                self.oldVol = self.StandVol
                self.MAI = self.StandVol / self.StandAge if self.StandAge > 0 else 0

                self.cRADint += self.RADint
                self.aRADint += self.RADint
                self.cGPP += self.GPP
                self.aGPP += self.GPP
                self.cNPP += self.NPP
                self.aNPP += self.NPP
                self.cCVI += self.CVI
                self.cLitter += self.lossWF
                self.cStemDM += self.incrWS
                self.aStemDM += self.incrWS
                self.cRainInt += self.RainIntcptn
                self.cTransp += self.Transp
                self.aTransp += self.Transp
                self.cEvapTransp += self.EvapTransp
                self.aEvapTransp += self.EvapTransp
                self.cRunOff += self.RunOff
                self.aRunOff += self.RunOff
                self.aSupIrrig += self.SupIrrig
                self.cSupIrrig += self.SupIrrig
                self.cLAI += self.LAI / self.monthsInStep

                # 输出频率控制
                current_time = f"{current_year}-{current_month:02d}"
                if self.output_frequency == "m":
                    # 月度输出：记录每个月的数据
                    results["Time"].append(current_time)
                    for var in output_vars:
                        if var in output_mapping:
                            results[var].append(output_mapping[var](self))
                        else:
                            print(f"Warning: Output variable '{var}' not supported.")
                            results[var].append(0)
                elif self.output_frequency == "a" and current_month == target_month:
                    # 年度输出：仅记录目标月份（种植月份）的数据
                    results["Time"].append(current_time)
                    for var in output_vars:
                        if var in output_mapping:
                            results[var].append(output_mapping[var](self))
                        else:
                            print(f"Warning: Output variable '{var}' not supported.")
                            results[var].append(0)

                current_month += 1
                if current_month > 12:
                    current_month = 1
                    current_year += 1

                if self.cRADint != 0:
                    self.cEpsilonStem = 100 * self.cStemDM / self.cRADint
                    self.cEpsilonGross = 100 * self.cGPP / self.cRADint
                if self.cEvapTransp != 0:
                    self.cWUE = 100 * self.cNPP / self.cEvapTransp
                if self.StandAge > self.ltAge0 + 0.99:
                    self.ltStemGR = (self.WS - self.ltWS0) / (self.StandAge - self.ltAge0)
                if self.opFrequency == 3 or monthCounter == 11:
                    if self.cLAI > self.LAIx:
                        self.LAIx = self.cLAI
                        self.ageLAIx = self.StandAge
                    if self.MAI > self.MAIx:
                        self.MAIx = self.MAI
                        self.ageMAIx = self.StandAge

        # 年度输出时，追加最终结果
        if self.output_frequency == "a":
            final_time = f"{self.planted_year + self.end_age}-{target_month:02d}"
            results["Time"].append(final_time)
            for var in output_vars:
                if var in output_mapping:
                    results[var].append(output_mapping[var](self))
                else:
                    print(f"Warning: Output variable '{var}' not supported.")
                    results[var].append(0)
        #logging.info(f"Step 2: Intermediate result {results}")
        return results

